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We present a detailed study of the ground-state magnetic structure of ultrathin Fe films on the surface of 
fee Ir(OOl). We use the spin-cluster expansion technique in combination with the relativistic disordered local 
moment scheme to obtain parameters of spin models and then determine the favored magnetic structure of the 
system by means of a mean field approach and atomistic spin dynamics simulations. For the case of a single 
monolayer of Fe we find that layer relaxations very strongly influence the ground-state spin configurations, 
whereas Dzyaloshinskii-Moriya (DM) interactions and biquadratic couplings also have remarkable effects. To 
characterize the latter effect we introduce and analyze spin coUinearity maps of the system. While for two 
monolayers of Fe we find a single-g spin spiral as ground state due to DM interactions, for the case of four 
monolayers the system shows a noncollinear spin structure with nonzero net magnetization. These findings 
are consistent with experimental measurements indicating ferromagnetic order in films of four monolayers and 
thicker. 

PACS numbers: 75.70.Ak, 71.15.Mb, 71.15.Rf 



I. INTRODUCTION 

By the end of the 20th century, in particular, since the birth 
of spintronics, thin films and nanostructures have gained in- 
creasing importance in industrial applications. The increas- 
ing demand for ultra-high density magnetic data storage de- 
vices had been one of the greatest driving forces of research 
and development involving nanostructures. With current hard 
disk technology approaching the superparamagnetic limit, the 
study of finite temperature magnetism in thin films and nanos- 
tructures is inevitable. Ab initio electronic structure methods 
give, in general, a very good description of the ground-state 
properties of solids. When trying to describe complex mag- 
netic structures or finite temperature magnetism these meth- 
ods are often used to generate parameters of spin models. 

It is by now widely accepted that relativistic corrections to 
the Heisenberg model, especially the Dzyaloshinskii-Moriya 
(DM) interaction, play an important role in determining the 
magnetic ground state of some systems.'"^ Moreover, higher 
order interaction terms (multiple spin interactions) also have 
to be considered in many cases^"^ to give an accurate descrip- 
tion of magnetism. Very recent studies indicate that such in- 
teractions may even lead to the formation of exotic states like 
magnetic skyrmion lattices.^ 

In the present work we demonstrate the use of spin mod- 
els from first principles for the description of magnetism, 
in particular, determining the ground-state spin configura- 
tion of thin film systems. We use the spin-cluster expan- 
sion (SCE) combined with the relativistic disordered local 
moment (RDLM) theory to obtain model parameters. The 
SCE-RDLM method has just been successfully applied to 
the IrMn3/Co(lll) interface, a prototype for an exchange 
bias system, to calculate exchange interactions and magnetic 
anisotropics.^ 



Here we go beyond the tensorial Heisenberg model by 
including biquadratic couplings in the spin Hamiltonian and 
present our results for thin films of Fe on the (001) surface 
of fee Ir. Previous theoretical studies have found that in case 
of a single Fe overlayer the favored magnetic configuration 
depends very strongly on layer relaxations, leading to the for- 
mation of spin spiral states near the experimental geometry. ' ' 
We study the effect of layer relaxations on the magnetic in- 
teractions for the case of a monolayer, as well as the effect 
of higher order interactions on the spin configuration. We 
also examine thicker films consisting of two and four layers 
of Fe to see whether the bulk ferromagnetism of Fe emerges 
with increasing film thickness. Experiments have shown that 
thin films of Fe consisting of 4 or more monolayers produce 
a ferromagnetic signal in magneto-optic KeiT effect (MOKE) 
measurements with an in-plane easy axis. The results of our 
simulations turn to be consistent with the main experimental 
findings. 



II. THEORY 

A. Spin model 

The most widely used ab initio methods to describe itin- 
erant magnetic systems rely on the adiabatic approximation 
separating fast single-electron spin fluctuations from the slow 
transversal motion of the spins. Furthermore, in the so- 
called rigid spin approximation it is assumed that longitudi- 
nal fluctuations of the local moments are negligible, so that 
the system of moments is characterized by a set of unit vec- 
tors {e} — {ei,..., e^} describing the orientation of each local 
moment. Then the grand potential i2({e}) of the system may 
be thought of as a classical spin Hamiltonian,''* which can be 
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used in numerical simulations to study the various magnetic 
properties of the system. 

For practical applicability, the energy must be parametrized 
in a simple yet meaningful way. The most common approxi- 
mation is a form of a generalized Heisenberg model, 



N ^ N 

i=l ^ ij=l 



(1) 



where K,- and J,y are the second order on-site anisotropy ma- 
trices and tensorial exchange couplings, respectively. The lat- 
ter can be meaningfully decomposed into an isotropic com- 
ponent jjj ~ TrJ,y73, an antisymmetric component jfj and 
a traceless symmetric part jfj. The isotropic component de- 
scribes a Heisenberg interaction, the antisymmetric part cor- 
responds to the DM interaction'^"'^ in the form of 



'iJfjej =Dij{eiXej) 



(2) 



and the final component contributes to the so-called two-site 
anisotropy. In this paper we go beyond the second order ex- 
pansion of Eq. (1) and include isotropic biquadratic interac- 
tion terms of the form — B,y(e, -ej)^. 

Even though a parametrization of the energy with a spin 
Hamiltonian possesses a much less direct connection to the 
magnetic properties of the system, it can be used well to 
simulate the magnetic behavior One method for obtaining 
the parameters of the spin Hamiltonian directly from first 
principles calculations is the so-called spin-cluster expansion 
(SCE) combined with the relativistic disordered local moment 
scheme (RDLM). 



B. Spin-cluster expansion 

The spin-cluster expansion developed by Drautz and 
Fahnle'^'"^ gives a systematic parametrization of the adiabatic 
energy surface. Up to two-spin interactions, the grand poten- 
tial may be expanded using real spherical harmonics as 

+ il I L 4fYdet)YL'iej), (3) 
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where the summations do not include the constant spherical 
harmonic function of the composite index {£,m) = (0,0). The 
coefficients in Eq. (3) are defined as 
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where vectors in lower index indicate restricted averages, i.e., 
uniform directional averaging has to be carried out with re- 
spect to every spin in the system not noted in the lower index. 



The terms of the spin Hamiltonian can be directly related 
to the terms of the SCE, for instance the isotropic biquadratic 
couplings can be expressed as 



Bij = - 



2 
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Clearly the key quantities of the SCE are the restricted direc- 
tional averages of the grand potential. 



C. Relativistic disordered local moment theory 

The DLM scheme gives a description of a magnetic system 
in accordance with the adiabatic approximation. Its imple- 
mentation within the Korringa-Kohn-Rostoker (KKR) theory 
was given by Gyorffy et a/.,''* with a relativistic generaliza- 
tion by Staunton et a/. '^'^^ Combining it with the SCE pro- 
vides a highly effective tool for determining the parameters of 
spin models. For a detailed presentation of the SCE-RDLM 
method see Ref 9. In the following we will review the most 
important features of the theory. 

The electronic charge and magnetization densities are de- 
termined from a self-consistent-field KKR calculation. In 
good moment systems the magnitude of local moments may 
be considered as independent from their orientation. For a 
given set of self-consistent potentials, charge and local mo- 
ment magnitudes, the orientations {?} of the local moments 
are accounted for by the similarity transformation of the 
single-site f-matrices. 



(8) 



where f,(ez) = is the f-matrix for a given energy, e 

(not labeled explicitly), with exchange field along the z axis, 
and R{ei) is the representation of the S0{3) rotation that trans- 
forms into e,. Underlines denote matrices in the {K,jj.) an- 
gular momentum representation. 

The coherent potential approximation (CPA) is employed 
to describe the magnetically disordered system. The strat- 
egy of the CPA is to substitute the disordered system with an 
effective (coherent) medium which is independent from the 
orientation of local moments, such that the scattering of an 
electron in the effective medium should resemble the average 
scattering in the disordered physical system. The scattering 
path operator of the effective medium is defined as 



-1 



(9) 



where double underlines denote matrices in site-angular mo- 
mentum space, is the matrix of structure constants, and 
is site diagonal. Using the excess scattering matrices defined 
as 
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(10) 



the single-site CPA condition can be formulated as 

Jd^eiX,{ei)=0. (11) 
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Within vaUdity of the magnetic force theorem the grand po- 
tential of the system can be expressed in terms of the excess 
scattering matrices and the related impurity matrices 



(12) 



for a given spin configuration and at zero temperature as 
= a. - - ^Im^ de lndetD,.(e,) 
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(13) 



where the constant term reads as 



deA^o(e) - - Im / de Indet t (e) , (14) 

7t 



No{e) being the integrated density of states of free electrons. 
Eq. (13) can be used to calculate restricted averages of the 
grand potential for the SCE, Eqs. (4)-(6). In particular the 
two-site expansion terms are expressed as 



J^f = --lmJ deJJ A'ci d'ej YL{e,)YL:{ej) 

X Trln [[-X,{ei)l.,jXj{ej)l,j^] , (15) 

implying that higher order two-site terms, such as biquadratic 
couplings, see Eq. (7), can be calculated just as easily as ten- 
sorial Heisenberg interactions. 
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III. RESULTS 
A. Computational details 

We used the screened Korringa-Kohn-Rostoker (SKKR) 
method' ''^^ to perform self-consistent calculations for fee 
bulk Ir and the layered Fe/lr(001) systems. The in-plane lat- 
tice constant for the fee lattice was chosen 2.715 A.'^ We used 
the local spin-density approximation parametrized according 
to Ceperley and Alder,^^ and we employed the atomic sphere 
approximation with an angular momentum cutoff of ^max — 3. 
A scalar-relativistic DLM description was used, correspond- 
ing to the paramagnetic state. Matrices of the effective CPA 
medium were determined to a relative error of 10"^ We used 
12 energy points on a semicircular path on the upper complex 
half -plane for the energy integrations, and 78 points were sam- 
pled in the irreducible wedge of the 2D Brillouin zone (BZ) 
for A;-integrations. 

The spherical integrations needed for the SCE, as in 
Eq. (15), were performed according to the Lebedev-Laikov 
scheme.""* The logarithm in Eq. (15) was expanded into a 
power series to avoid the phase problem due to the energy 
integration of the complex logarithm. 
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FIG. 1. (Color online) Isotropic couplings for various layer relax- 
ations in Fei/Ir(001) as a function of interatomic distance in units of 
the in-plane lattice constant, fl2d- 



B. Fei/Ir(001) 

Firstly we performed calculations for the case of an Fe 
monolayer. Along with the unphysical case of an unrelaxed 
Fe monolayer, systems with different layer relaxations were 
considered, namely 5%, 10%, 15% inward relaxations as well 
as the experimental value of -12% relaxation.'^ According 
to Kudrnovsky et a/.," we may expect a very strong depen- 
dence of exchange parameters on the layer relaxation, with 
a crossover of dominant Heisenberg couplings from ferro- 
magnetic (FM) to antiferromagnetic (AFM). The obtained 
isotropic couplings are shown in Fig. 1. 

While for the unrelaxed geometry all significant couplings 
are FM, with increasing inward relaxation the couplings for 
the three nearest neighbor (NN) shells become gradually 
AFM. For the case of experimental layer relaxation, the mag- 
nitude of the Heisenberg interactions is small and comparable 
up to four shells. While the tendency of the obtained curves 
is similar to those calculated by Kudrnovsky et al.,^^ there 
are remarkable differences. In particular, at the experimen- 
tal layer relaxation our largest isotropic couplings are 1 mRy 
smaller than those obtained in Ref. 11, corresponding to a 
difference by a factor of three. Since the TB-LMTO-SGF 
method is known to give results that agree very well with KKR 
we use, we performed a series of check calculations to verify 
the reliability of our results. Increasing the numerical preci- 
sion of the spherical integrations of the SCE, the convergence 
parameters of the CPA and the number of k points for the BZ 
integration all resulted in identical isotropic couplings within 
linewidth. Increasing the number of energy points from 12 to 
16 used for the energy integrations in the self-consistent cal- 
culations also leads to relative differences less then 2%. 

We also performed calculations with spin-orbit coupling 
scaled to zero^^ as well as with the parametrization of the 
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TABLE I. Calculated DM interactions and isotropic biquadratic cou- 
plings, see Eq. (7), in Fei/Ir(001). All values are given in mRy units. 



Layer relaxation 


0% 


-5% 


-10% 


-12% 


-15% 


1st shell 


0.023 


0.142 


0.244 


0.280 


0.331 


1 2nd shell 


0.156 


0.207 


0.209 


0.198 


0.170 


1st shell 


0.110 


0.109 


0.088 


0.077 


0.059 


ij 

2nd shell 


0.008 


0.012 


0.013 


0.013 


0.012 



LDA according to Vosko, Wilk and Nusair. The differences 
in the dominating exchange interactions from both sets of cal- 
culations were less then 0. 1 mRy indicating that such factors 
are insufficient to explain the quantitative difference between 
our interactions and those of Kudrnovsky et al. ' ' 

Due to the frustrated nature and small magnitude of the 
isotropic exchange interactions it is possible that DM inter- 
actions or biquadratic couplings affect the ground-state spin 
configuration. The magnitude of these interactions for the 
first two NN shells is shown in Table I. For higher relaxation 
values, in particular for the experimental one, the correction 
terms are indeed comparable to the isotropic couplings. Inter- 
estingly, in our case the biquadratic couplings all have posi- 
tive sign indicating that these interactions favor collinear spin 
configurations. This feature implies a competition between 
the frustrated AFM exchange couplings, the DM and the bi- 
quadratic interactions, possibly leading to complex ground- 
state spin configurations. 

The Fourier transform of the tensorial coupling matrices 
was calculated to obtain the mean field estimate, as explained 
in the Appendix. In the monolayer case the Fourier transform 
is a 3 X 3 tensor field defined on the two-dimensional BZ. The 
largest eigenvalue J{q) can be visualized as a surface in recip- 
rocal space. The mean field estimate states that the maximum 
points of the surface give the characteristic spatial modulation 
of the ordered magnetic state, to which the paramagnetic state 
is unstable. 

Reassuringly, the J{q) surface for the unrelaxed system has 
a pronounced maximum in the F point of the BZ, indicating 
FM ground state. As the layer relaxation is increased, the cor- 
responding surfaces are gradually depressed at the F point, 
giving way to new maxima in general points of the Brillouin 
zone for -10% and larger inward relaxations. The contour plot 
of the J{q) surface for experimental layer relaxation is shown 
in Fig. 2. The maxima are around q = (0.5,0.7)^ and sym- 
metry related points. This estimate suggests a complex non- 
collinear spin structure due to the frustrated and competing 
magnetic interactions detailed earlier. Note that there is an 
approximate degeneracy along lines connecting pairs of max- 
ima. The difference between the maximum values of the sur- 
face and its value at, for instance, q — (0.6,0.6)— is less then 
0.03 mRy. This implies the possibility that the actual ground 
state of the system is somewhere along these degenerate lines, 
but not necessarily in the precise maxima of the surface. 

The effect of relativistic tensorial interactions may be as- 
sessed by calculating the J{q) surface using only the isotropic 
part of the exchange tensors. In doing so, the shape of the 
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FIG. 2. (Color online) J(q) surface for Fei/Ir(001) with experimental 
relaxation. The color coding in units of mRy is shown on the right. 

resulting surface is the same, however its maxima are shifted 
closer to the zone boundary, to ^ = (0.45,0.85) ^ and related 
points. The degeneracy of the line connecting the maxima is 
remarkably decreased. Based on the mean field estimate, the 
ground-state spin configuration of Fei/Ir(001) is a spin spiral, 
the periodicity of which is largely affected by the DM interac- 
tions. 

We performed zero temperature Landau-Lifshitz-Gilbert 
spin dynamics simulations to verify the mean field prediction 
and to reveal the influence of biquadratic couplings on this 
noncoUinear spin structure. Due to the expected noncollinear 
structures we used a 128 x 128 lattice of spins along with free 
boundary conditions to prevent the appearance of spurious pe- 
riodicity. Interactions were included up to 7 two-dimensional 
lattice constants (fl2d)- For every relaxation two kinds of sim- 
ulations were performed, one with and the other without bi- 
quadratic coupling terms. Every simulation was initialized 
from a random spin configuration, the same for simulations 
with and without biquadratic terms. All simulations were con- 
tinued until there was no measurable difference (i.e., less than 
10^^ mRy) in the total interaction energy of the system. It 
should be noted that on-site anisotropy was not included in 
the simulations due to its small size (0.06 mRy in case of ex- 
perimental layer relaxation), however two-site anisotropy was 
included as the symmetric part of the J_.^ tensors. 

In the case of unrelaxed geometry, the simulations led to 
(out-of-plane) FM order as expected. As layer relaxation is in- 
troduced, the transition predicted by the mean field estimates 
is reflected in the spin dynamics simulations as well, as the 
ground states become complex spin-spirals for relaxations of - 
10% and above. The Fourier components of the emergent spin 
configurations agree very well with the mean field estimates. 
Interestingly, for -15% relaxation the spin configuration of the 
simulated system contains a huge number of domains with 
various orientations for the spin spirals, which is in very good 




FIG. 3. (Color online) Ground-state spin configurations of the 
Fe monolayer with experimental layer relaxation (a) without and 
(b) with biquadratic couplings. 



agreement with the degeneracy of the J{q) surface along the 
corresponding ^points. 

It is expected that biquadratic couplings can significantly 
affect the ground-state spin configuration when the tensorial 
Heisenberg interactions alone lead to the formation of non- 
collinear spin structures. This is the case for the experimental 
layer relaxation for which the resulting spin-configurations of 
the two simulations are shown in Fig. 3. The lattice Fourier 
transform of the spin configuration in Fig. 3a has a peak at 
q sa (0.719,0.495) — , which perfectly coincides with the nu- 
merical maximum of the corresponding J{q) surface. Interest- 
ingly the peak for the simulation including biquadratic terms 
is at ^ w (0.735,0.500) which is only very slightly differ- 
ent from the previous value. Even though the two spin struc- 
tures in Fig. 3 seem very different, their spatial modulation 
is almost identical. The main reason for this is that while in 
Fig. 3b there is a clear two row-by-two row periodicity along 
one direction, in Fig. 3a we can see a helical spiral along the 
same direction with 90° angle between adjacent spins. Both 
orderings have a period of 4fl2d- 

We may try to grasp the effect of biquadratic interactions 
based on the real space spin structure. From the form of the bi- 
quadratic coupling between two specific spins, — B,; (e,- • e^)^, 
it is clear that a coupling with positive sign favors collinear 
alignment, either parallel or antiparallel. It is motivated that 
positive biquadratic couplings acting on a noncollinear spin 
configuration will try to increase the coUinearity within the 
system. 




FIG. 4. (Color online) CoUinearity map for the approximate ground 
state of Fei/Ir(001) with experimental layer relaxation, (a) without 
and (b) with biquadratic couplings. The color coding in units of mRy 
is shown on the right. 



Let US define the coUinearity of two spins as 



co\l{e i,ej) 



arccos (e,- -ej) - 



(16) 



which is simply the deviation of the angle between the two 
spins from a right angle, normalized between and 1, 1 in- 
dicating most collinear aiTangement (i.e., where e,- • ej = ±1). 
Using this definition we may create a coUinearity map of any 
spin configuration, plotting at each site an averaged coUinear- 
ity defined as 



coll(;) = -y coll(e/,e;), 

^' J 

i'J) 



(17) 
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where the summation includes every first NN of site / and cor- 
respondingly Zi stands for the coordination number of site /. 
The reason for the consideration of only first NN sites is, be- 
sides practicality, the rapid spatial decay of biquadratic cou- 
plings, suggesting that this interaction is mainly sensitive to 
the first NN collinearity. 

Using the above definitions, the calculated collinearity 
maps from the two kinds of simulations are shown for ex- 
perimental relaxation in Fig. 4. The difference between the 
two maps is remarkable. At first glance we can see that the 
two configurations have a completely different domain struc- 
ture. While in Fig. 4a there are smoothly connected domains 
of low collinearity (0.22 on average), in Fig. 4b we can see 
neatly separated, homogeneous domains of high collinearity 
(0.67 on average). This is a clear indication that biquadratic 
couplings have a serious effect on the magnetic structure in 
this specific system. 
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C. Fe2/Ir(001) 

Martin et aO^ found using MOKE measurements that Fe 
films of 4 monolayers or thicker produce an in-plane FM sig- 
nal at room temperature. For the case of an Fe monolayer with 
reasonable layer relaxations the theoretically predicted spin 
spiral states have zero net magnetization, thus they would not 
provide a MOKE signal. To see if the bulk-like behavior of Fe 
emerges for thicker films, we extended our studies to cases of 
two and four monolayers of Fe on Ir(OOl). 

The geometries of Fe„/Ir(001), both for n — 2 and n — A 
were assumed according to Ref. 12, with the minor simplifi- 
cation that no layer relaxation was applied to the Ir substrate. 
As shown in Fig. 5, the isotropic Heisenberg couplings we 
obtained show a strong FM coupling between the two lay- 
ers, while the intralayer couplings are smaller, and in case of 
the subsurface layer the dominant ones are AFM. It should 
be mentioned that the DM interactions are one order of mag- 
nitude smaller then the isotropic terms, and the biquadratic 
couplings are even smaller by one order. 

The J{q) surface for the bilayer has its maximum near the 
F point, but not exactly at the zone center. There is a de- 
generate circle at \q\ k, 0.07— where the surface is maximal, 

fe 111 

predicting a long wavelength spin spiral. Considering only the 
isotropic part of the exchange tensors the corresponding sur- 
face has its maximum in the F point. This indicates that the 
DM interaction imposes a noncoUinear spin structure on the 
otherwise FM system. 

Since the number of neighbors within the same radius 
rapidly increases with the number of Fe layers, the spin dy- 
namics simulations for the multilayers were performed with 
in-plane system sizes of 64 x 64 sites. Fig. 6 shows the re- 
sulting spin configuration of the surface Fe layer for the sim- 
ulation including biquadratic couplings. Note that due to the 
strong FM interlayer coupling the subsurface layer displays a 
very similar pattern. 

The Fourier transform of the spin configuration indicates 
that the propagation direction of the spin spirals is not exactly 
diagonal, since the wave vector is q — (0.13,0.16)^ (and 



FIG. 5. (Color online) Isotropic Heisenberg couplings in Fe2/Ir(001) 
as a function of interatomic distance in units of the in-plane lattice 
constant, a2d- 
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FIG. 6. (Color online) Approximate ground-state spin configuration 
of Fe2/Ir(001), surface layer shown from the top. Arrows are colored 
according to z component, ranging from red (up) to blue (down). 



symmetry related points). The fact that q is not along any 
high symmetry lines is consistent with the degeneracy seen 
on the y(^) surface, however, the magnitude of the wave vec- 
tor is twice as large as that of the mean field estimate. This 
quantitative difference can mainly be attributed to the flaws 
of the mean field approach, but we can not rule out finite size 
effects in the spin dynamics simulation. 

The simulation without biquadratic interactions resulted in 
a very similar pattern indicating that these interactions are ir- 
relevant in thicker films due to increased isotropic Heisenberg 
couplings. Simulations performed using only the isotropic 
part of the Heisenberg tensors led to the appearance of large 
FM domains without the spin spiral pattern seen in Fig. 6. 
This result is corroborated by several previous observations 
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that the DM interaction may significantly affect the ground- 
state spin configuration, and even lead to the formation of spin 
spirals in thin film systems.'^'^ It is also interesting to note that 
a very similar labyrinth pattern was found in a monolayer of 
Mn on W(001).2 Composed of single-g' cycloidal spin spirals, 
the configuration shown in Fig. 6 does not have a net magnetic 
moment. 



D. Fe4/Ir(001) 

In the following, for the case of four monolayers of Fe, the 
Fe layers are going to be indexed from the substrate to the 
surface, so that layer 1 is closest to the Ir substrate and layer 4 
is the surface layer. Summarizing the calculated interactions, 
the interlayer isotropic Heisenberg couplings are strong and 
FM, with increasing magnitude towards layers closer to the 
surface. Interestingly the couplings between intralayer first 
nearest neighbors are mostly weakly antiferromagnetic, only 
the two layers closest to the substrate show couplings of con- 
siderable magnitude, above 1 mRy in absolute value. The cor- 
rections to the Heisenberg spin-model are weak. Only the first 
NN intralayer DM vectors in layer 1 are larger than 0.1 mRy, 
and the strongest biquadratic couplings are around 0.07 mRy. 
Based on these features it is probable that the magnetic behav- 
ior is dominated by isotropic couplings. 

The J{q) surface for Fe4/Ir(001) shown in Fig. 7 has its 
maximum in the F point, anticipating FM arrangement. If the 
ground state is indeed FM, it is worth looking at the mean 
field estimate for the Curie temperature. The maximum of the 
surface with 10.89 mRy corresponds to a mean field estimate 
of Tc — 573 K. As mean field approximations overestimate the 
stability of the ordered phases due to neglected fluctuations, 
the actual critical temperature is surely lower than this value. 
Still, it is possible that at room temperature the system is in 
the ordered phase. 

The spin dynamics simulations revealed that the ground 
state is more complicated than what the mean field approach 
predicts. The simulated system converged into a complex spin 
structure regardless of the presence of biquadratic coupling 
terms. In each layer the spins formed a cycloidal spiral of 
q = (0.41,0.41)^ rotating in a vertical plane, but with an ad- 
ditional in-plane FM modulation leading to nonzero net mag- 
netization, (e,) 7^ 0. While the wave vector of the spiral is the 
same within every layer due to the strong interlayer coupling, 
the magnitude of the FM Fourier component increases with 
distance from the substrate. For a quantitative comparison of 
the net magnetizations, the layer-averaged magnetizations cal- 
culated from a homogeneous domain containing 1178 spins 
are presented in Table II. The average magnetization, (e,), is 
clearly in-plane and monotonically increasing in magnitude 
towards the surface where it reaches a value larger than 0.7. 
Considering the actual size of each magnetic moment, also 
shown in Table II, the total average magnetic moment of the 
system is 1.27/iB- 

In contrast to bulk bcc Fe, the net magnetization points 
along the (110) direction, i.e., towards intralayer second- 
nearest neighbors, but it is indeed in-plane. Possibly in even 




FIG. 7. (Color online) J{q) surface for Fe4/Ir(001). The color coding 
in units of mRy is shown on the right. 



TABLE 11. Layer-resolved averaged magnetizations, (e,), and atomic 
magnetic moments in a given domain of Fe4/Ir(001) containing 1 178 
spins. 







magnetic moment [jls] 


layer 1 


(0.26,0.26,0.01) 


1.98 


layer 2 


(0.37,0.37,0.01) 


2.02 


layer 3 


(0.51,0.53,0.00) 


1.60 


layer 4 


(0.54,0.57,0.00) 


2.71 



thicker films we would see FM order with easy axis along the 
(100) direction, as was found experimentally. 



IV. CONCLUSIONS 

We used the recently developed SCE-RDLM method to ob- 
tain spin Hamiltonians from ab initio calculations, going be- 
yond the anisotropic Heisenberg model by including isotropic 
biquadratic interaction terms. The obtained interaction pa- 
rameters allow for a detailed investigation of the magnetic 
ground state via a mean field approach and atomistic spin dy- 
namics simulations. We presented results for Fe thin films 
on the (001) surface of a semi-infinite Ir substrate. For the 
case of an Fe monolayer we found that layer relaxations dras- 
tically rearrange the interaction landscape, leading to the ap- 
pearance of complex non-coUinear spin structures. Relativis- 
tic corrections, in particular, Dzyaloshinskii-Moriya interac- 
tions, are needed to be taken into account as anisotropic cou- 
plings largely affect the ground state. Spin dynamics simu- 
lations also revealed that including biquadratic interactions to 
the spin model significantly alters the ground state spin con- 
figuration by favoring a more collinear state. This serves as 
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an instructive warning that in some systems the usual Heisen- 
berg model might give an insufficient description of magnetic 
properties. 

For thicker films of two and four monolayers we found 
that biquadratic couplings are irrelevant in determining the 
ground-state spin configuration due to the large magnitude of 
the Heisenberg interaction. While the bilayer system still pro- 
duces a single-^ spin spiral as ground state due to the DM in- 
teraction, in the quadrilayer system there is nonzero net mag- 
netization superimposed on a cycloidal spin structure. This 
finding may be consistent with experimental results showing 
a ferromagnetic signal in MOKE measurements above four 
monolayers of Fe deposited on Ir(OOl). 
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Appendix: Mean field paramagnetic spin susceptibility 

For a system of spins {e} governed by the grand poten- 
tial £l({e}), the best variational mean field trial Hamiltonian 
Q.Q {{e}) = — Y^ihi^i, i.e., the one that minimizes the free en- 
ergy is given by 

hi = -^Jd^eiei{ao)l, (A.l) 

where the angle brackets now denote thermodynamic aver- 
aging with respect to the mean field probability distribution 
Po({e,}) — exp [— /3i2o({e,})] /Zq, Zq being the corresponding 
canonical partition function. 

For our spin model extended with an inhomogeneous exter- 
nal field, 

N 2 w 

i=l ij=l 

I%(?r?,)'-f /ir-?. (A.2) 

i,j=l i=l 

i¥j) 

the Weiss field reads as 

hi = hf^ + h^''\ (A.3) 

where the molecular field induced by the interaction takes on 
the familiar mean field form. 



rrij = (cj) being the average magnetization at site j. Note 
that neither the on-site anisotropy nor the biquadratic cou- 
plings contribute to the Weiss field. 

The spin susceptibility, defined as j"^ = dmf /dh^^^'^ may 
be related to the susceptibility of the non-interacting spin sys- 
tem X^j"^ — dmf/dh^ as 

X<7 = Xij + L xUklXlj- (A.5) 

For layered systems with only two-dimensional translation in- 
variance we may separate site indices according to a layer in- 
dex and an in-plane index, 

XiJJj^Xu.ij+ L X^.ik^KL.kiXujj, (A.6) 

K.L.kJ 

(K,kmu) 

where capital indices denote layers. Due to in-plane transla- 
tion invariance we may introduce the two-dimensional Fourier 
transform of the quantities, for instance 

X/y(^)=Exw.'0e-'5^', (A.7) 
Ri 

which might be thought of as blocks of matrices in layer and 
Cartesian space, 

[xm,j = Xij{q)- (A.8) 

Using this notation we may easily invert Eq. (A.6) and ar- 
rive at 

X{q) = (l - x\q)m) X\q) , (A.9) 

indicating an enhancement of the spin susceptibility due to the 
Heisenberg interaction. In the paramagnetic limit this simpli- 
fies to 

x{q)={^kBTl-mY\ (A.IO) 

Annealing the system from the paramagnetic phase there will 
be a temperature, T^^a, where x{q) becomes singular, indicat- 
ing that the paramagnetic phase is unstable to the formation 
of some magnetically ordered state at that temperature. This 
transition temperature Toid is given by the condition 

Toid = max I eigenvalues of J(^) | 

= -^max/(^), (A.ll) 
iks q 

and the corresponding q where this happens gives the charac- 
teristic wave vector of the ordered phase (where the q vectors 
are sought for in the two-dimensional BZ). In Eq. (A.ll) we 
introduced the notation J{q) for the maximal eigenvalue of 
J(^) at a given wave vector. We may gain valuable insight re- 
garding magnetic ordering in the system by plotting this scalar 
function against the points of the BZ. 



9 



* deak.andris@gmail.com 

^ M. Bode, M. Heide, K. von Bergmann, P. Ferriani, S. Heinze, 
G. Bihlmayer, A. Kubetzka, O. Pietzsch, S. Bliigel, and 
R. Wiesendanger, Nature 447, 190 (2007). 

^ P. Ferriani, K. von Bergmann, E. Y. Vedmedenko, S. Heinze, 
M. Bode, M. Heide, G. Bihlmayer, S. Bliigel, and R. Wiesen- 
danger, Phys. Rev. Lett. 101, 027201 (2008). 

^ L. Udvardi, A. Antal, L. Szunyogh, A. Buruzs, and P. Weinberger, 
Physica B 403, 402 (2008). 

A. Antal, B. Lazarovits, L. Balogh, L. Udvardi, and L. Szunyogh, 
Phil. Mag. B 88, 2715 (2008). 
^ A. Antal, L. Udvardi, B. Ujfalussy, B. Lazarovits, L. Szunyogh, 
and P Weinberger, J. Magn. Magn. Mater. 316, 118 (2007). 

* B. Hardrat, A. Al-Zubi, P. Ferriani, S. Bliigel, G. Bihlmayer, and 
S. Heinze, Phys. Rev. B 79, 094411 (2009). 

^ S. Lounis and P H. Dederichs, Phys. Rev. B 82, 180404(R) 
(2010). 

^ S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, 
R. Wiesendanger, G. Bihlmayer, and S. Bliigel, Nature Physics 7, 
713 (2011). 

^ L. Szunyogh, L. Udvardi, J. Jackson, U. Nowak, and R. Chantrell, 

Phys. Rev. B 83, 024401 (2011). 
" L. Udvardi, L. Szunyogh, K. Palotas, and P. Weinberger, Phys. 

Rev. B 68, 104436 (2003). 
^ J. Kudmovsky, F. Maca, I. Turek, and J. Redinger, Phys. Rev. B 

80, 064405 (2009). 



V. Martin, W. Meyer, C. Giovanardi, L. Hammer, K. Heinz, 
Z. Tian, D. Sander, and J. Kirschner, Phys. Rev. B 76, 205418 
(2007). 

'■^ V. P. Antropov, M. L Katsnelson, B. N. Harmon, M. van Schilf- 
gaarde, and D. Kusnezov, Phys. Rev. B 54, 1019 (1996). 
B. L. Gyorffy, A. J. Pindor, J. B. Staunton, G. M. Stocks, and 
H. Winter, J. Phys. F: Met. Phys. 15, 1337 (1985). 

15 L Dzyaloshinskii, J. Phys. Chem. Solids 4, 241 (1958). 
T. Moriya, Phys. Rev. 120, 91 (1960). 

1^ R. Drautz and M. Fahnle, Phys. Rev. B 69, 104404 (2004). 
R. Drautz and M. Fahnle, Phys. Rev. B 72, 212405 (2005). 

1^ J. B. Staunton, S. Ostanin, S. S. A. Razee, B. L. Gyorffy, L. Szun- 
yogh, B. Ginatempo, and E. Bruno, Phys. Rev. Lett. 93, 257204 
(2004). 

J. B. Staunton, L. Szunyogh, A. Buruzs, B. L. Gyorffy, S. Ostanin, 
and L. Udvardi, Phys. Rev. B 74, 144411 (2006). 
L. Szunyogh, B. Ujfalussy, P. Weinberger, and J. Kollar, Phys. 
Rev. B 49, 2721 (1994). 

R. Zeller, P. H. Dederichs, B. Ujfalussy, L. Szunyogh, and 

P Weinberger, Phys. Rev. B 52, 8807 (1995). 
23 D. M. Ceperley, and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 
2'* V.L Lebedev, and D.N. Laikov, Doklady Mathematics 59, 477 

(1999). 

25 H. Ebert, H. Freyer, A. Vernes, and G.-Y. Guo, Phys. Rev. B 53, 
7721 (1996). 

2*' S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 
(1980). 



